clear
cd "C:\Users\JGW99\Dropbox\Forest\Data"

***************************************************************
*Creating biomass 5-county dataset from Yucheng
**************************************************************
/*
use temp\Biomass2011_NDVI_Bedford.dta
rename name taxidnum
gen fips=42009

append using temp\Biomass2011_NDVI_Centre.dta
replace fips=42027 if fips==.

append using temp\Biomass2011_NDVI_Huntingdon.dta
replace fips=42061 if fips==.

append using temp\Biomass2011_NDVI_Somerset.dta
replace fips=42111 if fips==.
replace taxidnum = mapnum if taxidnum==""

append using temp\Biomass2011_NDVI_Potter.dta
replace fips=42105 if fips==.
replace taxidnum = map_number if taxidnum==""

keep taxidnum fips v_biomassjks_rf_predicts_forest v_canopyheight_mean_forest v_treecover_percent_forest 

*Dropping repeated taxidnum within a county
duplicates drop taxidnum fips, force

save temp\Biomass2011Only, replace

*Bringing in 2007 NDVI dataset
gen year=2007
qui: merge 1:1 taxidnum fips year using temp/NDVI_landtr_5.dta, keep(3) keepusing(ndvi_* year) nogen 

*Bringing in forest acres variable
merge 1:1 taxidnum fips using ForestAcres, keep(3)

*Keeping only 1+ acre forested parcels
keep if forestAcres>=1

*renaming
rename v_canopyheight_mean_forest canopy_height_f
rename v_treecover_percent_forest treecover_forest
rename ndvi_landtr_mean_forest ndvi_forest

save temp\Biomass2011, replace
*/
*******************************
*Analysis
*******************************

use temp\Biomass2011

*Scatterplots

*Biomass minimum
drop if biomass_acre_f==0

gen biomass_ha=biomass_acre*2.47

*Estimation

gen ndvi_forest2=ndvi_forest^2 
gen ndvi_forest3=ndvi_forest^3 

/*Takeaway from non-linear estimation is that the slope is positive over all relevant ranges of NDVI*/

correlate biomass_ha ndvi_forest if ndvi_forest>800 & biomass_acre_f<95

regress biomass_ha ndvi_forest ndvi_forest2 if ndvi_forest>800 & biomass_acre_f<95, robust
xi: regress biomass_ha ndvi_forest i.fips if ndvi_forest>800 & biomass_acre_f<95, robust





